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ABSTRACT 


A method  was  developed  for  computing  the  displacement  spec- 
trum of  earthquake  dislocations  and  presumed  explosions  with  greater  accur- 
acy. Measurements  of  corner  frequencies  and  displacement  amplitudes  are 
considerably  less  ambiguous  than  those  derived  from  conventional  spectra. 

The  relatively  narrowband  system  response  is  apparently  correctly  accounted 
for  in  the  equations  for  ground  displacement,  so  that  the  displacement  spectra 
are  in  principal  valid  up  to  10  Hz.  Presently,  roundoff  errors  compounded  by 
estimation  errors  due  to  the  coda  of  the  events  limit  the  accuracy  and  validity 
of  spectra  at  low  frequencies  less  than  approximately  1 Hz.  It  iB  anticipated 
that  with  further  development  of  the  method,  valid  displacement  spectra  can 
be  measured  at  lower  frequencies. 

Discriminants  between  earthquakes  and  explosions  were  derived 
from  measurements  of  corner  frequency  and  displacement  amplitudes  taken 
from  the  spectra  of  events.  These  discriminants  appeared  to  be  at  least  as 
effective  as  other  short-period  spectral  discriminants  such  as  spectral  split- 
ting, spectral  moments,  and  spectral  magnitudes.  The  discriminants  utilized 
measurements  of  spectra  which  are  accurate  only  at  frequencies  greater  than 
1 Hz.  Additional  discriminants  based  on  low  frequency  characteristics  of  the 
displacement  spectra  may  be  possible  in  the  future  when  the  accuracy  of  dis- 
placement spectral  calculation  is  improved. 

The  spectral  measurements  of  corner  frequencies  and  displace- 

1/2 

ment  amplitude  were  used  to  compute  magnitudes  assuming  E scaling, 
where  E is  the  radiated  seismic  energy.  Comparison  of  spectral  measure- 
ments of  magnitudes  at  regional  distances  to  teleseismic  time  domain  mea- 
surements of  earthquake  magnitudes  by  the  seismic  network  indicates  that 


such  measurements  could  be  made  with  unexpected  precision  (approximately 
0.  17  standard  deviations  of  magnitude).  The  results  from  teleseismic  mea- 
surements of  P-waves  were  less  precise  (approximately  0.29  standard  de- 
viations). For  earthquakes,  magnitudes  derived  from  spectral  parameters 
are  consistent  with  the  event  mfa  derived  from  time  domain  measurements. 
For  explosions,  magnitudes  derived  from  spectral  parameters  are  consis- 
tently higher  than  those  derived  from  time  domain  measurements  by  a con- 
stant amount.  This  constant  positive  difference  between  the  spectral  and 
time  domain  magnitude  measurement  indicates  events  with  greater  high  fre- 
quency radiated  energy.  This  magnitude  difference  along  with  measurement 
of  the  spectrum  corner  frequency  appears  to  be  the  best  means  of  discrimina- 
ting explosions  from  earthquakes  based  on  the  high  frequency  characteristics 
of  the  events. 


t U • f' 1 'V  Advanced  Research  Projects  Agency  nor  the  Air  Force 
Technical  Applications  Center  will  be  responsible  for  information  contained 
herein  which  has  been  supplied  by  other  organizations  or  contractors,  and  this 
ocument  is  subject  to  later  revision  as  may  be  necessary.  The  views  and  con- 
clusions presented  are  those  of  the  authors  and  should  not  be  interpreted  as 
necessarily  representing  the  official  policies,  either  expressed  or  implied,  of 
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SECTION  I 
INTRODUCTION 


The  purpose  of  this  study  is  to  further  investigate  a possibility 

of  deriving  more  effective  short-period  seismic  discriminants.  In  the  past. 

short-period  discriminants  have  been  used  to  confirm  discrimination  results 

of  more  reliable  methods  such  as  M versus  m and  focal  depth  estimates. 

s b 

An  effective  set  of  short-period  discriminants  would  be  useful  in  its  own  right. 
The  detection  capability  for  short  period  data  is  in  most  cases  better  than 
for  long-period  data.  Also,  the  cost  of  installation  and  operation  of  a world- 
wide system  for  acquisition  of  long-period  data  is  considerably  greater  than 
that  for  short-period  data.  Therefore,  it  is  worth  directing  research  toward 
establishing  the  limit  of  applying  short-period  data  to  discriminating  explosions 
from  earthquakes.  Knowledge  of  this  limit  will  assure  that  the  final  system 
configuration  implemented  will  provide  the  most  cost  effective  seismic 
monitoring  of  nuclear  explosions. 


SECTION  II 


THEORETICAL  BACKGROUND  AND  PREVIOUS  EXPERIMENTAL  RESULTS 


A. 


THEORETICAL  BACKGROUND 


Briefly  described  is  a theoretical  basis  for  applying  measure- 
ments of  corner  frequency  and  other  spectral  parameters  to  the  problem  of 
discrmvnating  between  explosions  and  earthquakes.  Brune  (1970)  presented 
a model  for  the  spectrum  of  shear  waves  from  a dislocation.  His  formula  for 
the  average  spectrum  is: 


<n(f)>  = <R0c,>2fF(e) 


„ 1/f 

~ , . c 


1+  (f/f  )‘ 

C 


(ii-i) 


where 

and 
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F (c)  = | ^2-2eJ£  l-cos  (1.21  e f/f  ) J + e 2 | 

° - stress 

(3  = shear  velocity 

l1  = rigidity 

fc  = corner  frequency 

c = fractional  stress  drop 

= higher  corner  frequency  due  to  effect  of  fractional 
stress  drop 

- average  of  radiation  pattern  of  shear  waves. 

Figure  II-l  shows  a spectrum  calculated  according  to  Equation  II-l. 

If  the  stress  drop  from  a dislocation  is  100%  complete,  e =1. 

Then  the  spectrum  is  flat  nearly  up  to  the  corner  frequency  (f  ) and  falls 
- 2 c 
off  as  f above  the  corner  frequency.  Fore«l,  the  spectrum  is  flat  up  the 

corner  frequency  above  which  it  falls  oL  as  f"1.  At  higher  frequencies  its 

slope  gradually  increases  until  at  frequencies  greater  than  f = f /e  the 

o c '* 


II-l 


-y 

spectrum  falls  off  very  sharply  as  f ° (2<y  <10),  as  a s-de  lobe  of  F(<)  is 

o ' 

approached. 


Hanks  and  Wyss  (19/2)  presented  the  theoretical  spectrum  of 

P v/aves  from  an  earthquake.  The  functional  form  of  their  result  is  similar 

to  that  of  Brure  although  the  definition  of  some  parameters  such  as  f and 

c 

R.  . )>  is  different. 

80 


The  expression  for  the  radiated  energy  from  earthquakes  is 
(Hanks  and  Wyss  1972). 


(II-2) 


by  a straight  forward  application  of  Parseval's  theorem  to  the  volume  inte- 
gration of  kinetic  energy  density.  Equation  II-2  is  valid  for  shear  or  com- 
pressional  waves  under  any  model.  For  Brune's  model,  given  an  arbitrary 
exponent  -y  on  the  frequency  above  the  corner  frequency  fc,  the  radiated 
energy  of  P or  S waves  is  given  by  choosing  the  appropriate  factors  in: 


E (P.s)  • (<r,/3 ) • 


2 2 nJ (P’S>  , , 

PR  (Zrr)2  {]  (P.  S)  (i  +- i~)  . 


<R0gi  (P,  S)  > 


2 y-3 


(II- 3 ) 


Here 

I(P)  = 4 7r / 1 5 . 1(S)  = 2 jr/15 

(a, (3)  = compressional  velocity  or  shear  velocity 

R = the  distance  in  kilometers 

^ = is  the  density. 

Py  equating  the  energy  derived  empirically  by  integrating  the 
spectra  to  that  derived  in  Equation  (II- 3)  it  is  possible  to  solve  for  an 


II- 3 


energy"  cornc-  frequency,  f^.  This  is  an  alternative  method  of  measuring 
corner  frequency  which  involves  not  only  the  corner  frequency  f but  also 
other  details  of  the  spectrum.  Assuming  y = 2 we  can  solve(lI- 3)  for  the  energy 
corner  frequency  for  P waves. 

dfl  . (II-4) 


I he  constant  . 03  was  obtained  using  3.  32  for  the  density  and  U 

i/ed  amplitude  spectral  density  taken  as  1.0  below  the  measured 
quency. 


is  the  normal- 
corner  fre- 


(1970)  is: 


Another  theoretical  scaling  law  found  by  Thatcher  and  Hanks 


Log  (2  r)  = 


|ml  + 2.9 


2 

3 AG 


where  Aa  is  the  stress  drop  and  r is  the  radius  of  the  dislocation. 

tionship  for  stress  drop,  where  12^  is  corrected  for  attenuation, 

the  corner  frequency  in  Drune's  model,  is: 


The  rela- 

and  f is 
c 


LogAa  = Log  n ' + 3 Log  f + constant. 

O c 

Assuming  that  m(j  and  scale  similarly  with  energy,  we  have 

2 2 

(2r)  = y mb  - - Log  ^ . 2 Log  f + constant 

where  mb  is  a network  magnitude.  By  taking  the  constant  equal  to 
(2 r ) is  proportional  to  the  source  dimension. 


01-5) 


zero,  Log 
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As  a resul*  of  the  response  of  the  seismic  system,  a layered 
earth  and  absorption,  the  conventional  magnitude  probably  reflects  the  seismic 
source  wavelet  at  frequencies  between  0.7  and  1.  5 Hz.  Thus  elevated  spectral 
magnitudes  ( or  radiated  seismic  energies)  reflect  the  combined  effect  of 
higher  corner  frequencies  and  high  frequency  peaks  in  the  spectrum  of  the 
signal.  The  presence  of  the  latter,  although  indicated  by  the  spectra  of  some 
events,  cannot  be  considered  well  established  due  to  the  present  accuracy  of 
conventional  spectral  computations  at  frequencies  less  than  the  order  of  1 Hz. 


B.  PREVIOUS  EXPERIMENTAL  RESULTS 


Using  a data  base  of  136  local  California  events  observed  by 

a local  network  of  seismometers,  Thatcher  and  Hanks  (1973)  computed  local 

magnitude  and  other  source  parameters  from  measurements  of  the  long-period 

spectral  displacement  amplitude  fl  and  the  corner  frequency  f of  shear  waves 

o c 

from  local  events.  Their  emperical  scaling  law  was  : 

M = log  13  + ~ log  f +7.2. 

L o 2 c 


This  result  combined  with  radiated  energy  computed  from  Equation  II - 4 led 
to  the  following  relationship  be  tween  radiated  energy  and  Mj^with  good  corr- 
espondence with  the  observed  data  : 


Log  E =2.0  M + 8.  0 . 

° s L 

The  local  magnitude  is  seen  to  be  proportional  to  log  E 
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Following  these  results,  it  has  been  assumed  here  that  m 's 

b 

for  regional  and  teleseismic  events  can  also  be  estimated  from  measurements 

of  J2  corrected  for  attenuation  with  distance  and  f , and  are  also  proportional 

° 1/2  c 

to  log  (Eg  ).  That  is: 


m,  = log  £2  + — log  f + constant  . 

bs  o 2 c 


(II  ~ 7) 


Analysis  of  data  by  Richter  (1958)  shows  that  this  assumption  is  close  if  not 
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exact.  The  validity  of  the  assumption  can  also  be  assessed  by  the  results 
obtained  here. 


SECTION  III 

SPECTRAL  ESTIMATION  TECHNIQUES 

A.  THEORETICAL  DEVELOPMENT 

Several  methods  can  be  used  to  find  the  earth  displacement 
spectral  density.  The  conventional  method  is  to  compute  the  power  as  a 
function  of  frequency  using  the  discrete  Fourier  transorm.  (DFT)  on  sampled 
data  and  to  smooth  the  result. 

There  can  be  problems  in  interpreting  the  results  of  spectra 
computed  this  way.  In  particular,  this  method  implictly  assumes  that  the  data 
outside  the  transform  interval  are  zero,  which  is  a less  reasonable  assump- 
tion for  lower  signal  to  noise  ratios.  Also,  due  to  truncation  of  the  signal, 
leakage  occurs  from  the  spectral  estimate  at  one  frequency  into  other  f e- 
quencies  of  lower  power.  These  effects  cause  large  errors  of  unknown 
magnitude.  For  seismic  data  these  problems  are  intensified  by  the  fading  ran- 
dom character  of  the  recorded  signal,  and  such  errors  are  magnified  at  low 
frequencies  when  corrections  are  applied  for  the  system  response. 

The  seismic  signal  can  be  conceived  of  as  a deterministic 
wavelet  convolved  with  a random  function  representing  multiple  delayed 
sources,  multiple  paths  through  the  medium,  and  reverberations  due  to 
scattering.  Even  direct  transmissions  arriving  after  the  first  can  change  in 
shape  due  to  effects  of  propagation.  Thus,  there  can  be  expected  large  vari- 
ence  from  the  assumption  of  time  stationarity  of  the  signal  waveform  and  the 
assumption  of  a stationary  random  process  for  coda  amplitudes.  To  obtain 
the  most  efficient  and  accurate  statistical  estimate  of  the  signal  which 
contains  source  effects  ilone  it  is  useful  to  estimate  the  average  wavelet 
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corresponding  to  the  widest  possible  range  of  emergence  angle  from  the 
source  and  paths  through  the  medium.  To  obtain  this  result,  the  maximum 
entropy  method  uses  the  empirical  autocorrelation  function  as  input  data. 

It  computes  from  a large  sample  of  data  the  smallest  possible  number  of 
points  needed  to  lcpresent  the  average  seismic  signal. 

The  maximum  entropy  method  (Barnard,  1969)  uses  a specified 
number  of  lags  of  the  autocorrelation  function  0 as  input  data,  and  computes 
the  best  prediction  error  filter  r in  a least  squares  sense. 


[*]  £ = E El=  b , 0,  0 . . . o]  . 


(HI-1) 


Here 


is  the  Toplit.6  autocorrelation  function  matrix 


M- 


0(0)  . . 0(1)  . . . 0(N) 

0(i)  0(0)  . . . tf(N-l) 

0 (N)  0(N  -1)  . . . .0(0) 


(III-2) 


and  Eq  is  the  energy  of  the  input  times  series.  The  filter  r predicts  the 
values  of  0 outside  the  data  interval  and  thus  avoids  the  problems  associat- 
ed with  assuming  that  they  are  zero  there.  Problems  arising  from  spectral 
windowing  are  reduced  to  the  extent  that  all  of  the  values  of  the  autocorrela- 
tion can  be  predicted.  Furthermore,  since  the  prediction  filter  is  based  on 
that  part  of  the  input  which  is  deterministic  and  hence  predictable,  it  can  be 
used  to  estimate  the  spectrum  of  the  seismic  wavelet,  as  follows  : 

The  z transform  of  J*  = [7.  , 7^  ..  . 7J  is  £ T z”  , 

n=  1 


and  the  absolute  value  of  the  spectrum  obtained  by  setting  Z°=  e*nAta) 

z transform  is  the  desired  amplitude  spectrum.  At  is  the  time  between  sam 
pled  data  points. 
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Another  more  accurate  method  of  solving  for  the  prediction 
error  filter  is  provided  by  ar  algorithm  refered  to  as  the  Burg  technique, 
and  is  described  by  Ulrych  (1973).  This  method  minimizes  the  effect  of 
errors  in  the  empirical  correlation  function  by  estimating  the  prediction 
error  filter  coeffecients  directly  from  the  data.  Comparisons  of  the  accuracy 
of  spectra  computed  conventionally  by  the  DFT  and  by  the  maximum  entropy 
method  are  given  by  King  et  al.  (1974). 

Another  formulation  of  the'  maximum  entropy  method  was  used 
here  to  compute  the  displacement  amplitude  spectrum.  It  can  be  shown  that 
the  above  formulation  for  the  prediction  error  filter  coeffecients  is  an  approx- 
imation even  if  its  assumptions  are  satisfied.  Thevalidity  of  its  results 
depend  on  the  amount  of  data  used  to  estimate  the  filter.  The  required 
record  length  to  achieve  specified  accuracy  is  related  to  the  distance  of  the 
closest  roots  of  the  z transform  of  the  prediction  error  filter  coeffecients 
from  the  unit  circle  in  the  z plane.  The  z transform  of  the  spectrum  under 
Brune's  model  can  contain  multiple  roots  just  inside  the  unit  circle.  It  is 
therefore  preferable  to  obtain  an  exact  solution  for  the  inverse  autocorre- 
lation function,  A,  rather  than  estimating  it  by  inverting  the  z transform  of 
the  prediction  error  filter.  This  is  done  by  solving 

Ma.Q  Q_t=  [o,0,0,,Eo...,0.0,  o]  ( III  - 3 ) 

where  A1  = a ...a  ...  a 1 and  a.(z)  = 7Tz)  7 (z.)  . 

L n o n J i l i 

It  can  be  shown  that  this  solution  ^A  is  exact  for  any  length  of  sampled  data. 

This  is  true  even  if  the  poles  of  the  z transform  of  the  prediction  error  filter 

lie  on  the  unit  circle.  To  find  the  spectral  amplitude,  the  cosine  transform 

of  the  vector  is  inverted  and  the  square  root  extracted.  Then  E is  found 

o 

by  convolving  the  normalized  inverse  auto<  ot  relation  A with  the  empirical 
autocorrelation  function.  By  extending  the  convolution  beyond  the  length  of 
the  input  time  series,  the  error  associated  with  the  prediction  of  the  auto- 
correlation function  at  these  lags  can  be  estimated.  This  error  is  expressed 
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as  a percentage  of  the  value  at  zero  lag,  E 

Given  the  energy  Eq  and  the  normalized  inverse  autocorrela- 
tion function  A the  amplitude  spectrum  is: 


n"  “ (mg-sec)  (HI-4) 

1+2  2 2*nfAt  a cos  2jnfr2 

n=l  n 1 

where  At  is  the  sample  time.  This  spectrum  pertains  to  that  part  of  the 

empirical  autocorrelation  function  which  can  be  predicted.  For  that  reason 

the  method  is  limited  by  the  presence  of  unpredictable  coda,  ambient  noise, 

and  in  some  cases  by  zeros  of  the  system  response  or  in  the  signal  wavelet, 

especially  those  which  lie  on  or  near  the  unit  circle. 


n(f)  = 


Just  such  a problem  arises  when  computing  the  displacement 
spectrum  of  seismic  signals  as  actually  recorded.  The  seismometer  response 
transform  can  be  closely  represented  by  ti  e ratio  of  polynomial  z transforms 
P(z)  / Q(Z).  The  poles  of  the  system  response  are  easily  removed  by  multiply- 
ing the  z transform  of  the  observed  signal  by  Q(z),  which  can  be  shown  to  be 
a stable  procedure  and  to  increase  the  signal  to  noise  ratio  by  enhancing  high 
frequencies.  On  the  other  hand  the  zeros  of  the  system  response  transform 
at  NORSAR  involve  a 3 rd  order  root  on  the  unit  circle.  Their  removal  b} 
synthetic  division  by  P(z)  is  an  unstable  operation  which  enhances  roundoff 
errors  and  greatly  reduces  the  signal  to  noise  ratio  by  enhancing  low  frequency 
energy.  Therefore,  the  equations  for  the  maximum  entropy  spectrum  are 
modified  to  include  the  zeros  of  the  system  response  as  follows.  A modified 
inverse  autocorrelation  function  B is  solved  by: 

m-i  yT*  ■„[•••*,.  v *v-] 

— " [v  bn-r",b0’”-b„-l’  bnJ 
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where  R * [.  . . R]t  RQf  Rj,.  . . J is  the  autocorrelation  function  of  the  sys- 
tem response  and  B 1 is  the  inverse  autocorrelation  function  of  the  earth  dis- 
placement corrected  for  instrument  response.  The  equation  for  the  displace- 
ment spectrum  is  analogous  to  (111-4). 


n<f)  = 


E1/2  At 
o 


n 


1 + 2 b cos  2 7rnfAt 


n=  1 


1/2 


(m^i-sec). 


(HI -5) 


To  be  as  objective  as  possible,  data  windows  to  compute  the 
autocorrelation  function  are  automatically  determined.  The  initial  point 
of  a signal  gate  is  chosen  by  visual  inspection  of  the  trace.  Two  windows 
beginning  there  are  automatically  detected. 

The  first,  called  the  signal  window,  is  picked  by  calculating 
the  power  in  a fixed  length  gate  which  moved  down  the  signal  from  the  initi.  1 
point.  The  point  at  which  this  power  decreased  for  some  preassigned  number 
of  points  was  chosen  as  the  end  of  the  signal  gate. 

The  coda  gate  was  chosen  by  calculating  the  signal -to-noise 
ratio  based  on  noise  preceeding  the  signal  and  a moving  signal  gate  as  before. 

When  this  signal -to -noise  ratio  dropped  below  2.0,  the  end  of  the  coda  was 
declared. 

Spectra  calculated  using  this  detector  are  relatively  in- 
sensitive to  its  parameters.  This  procedure  minimizes  the  influence  of 

varying  coda  lengths  in  the  determination  of  the  signal  pulse  displacement 
spectrum. 


B. 


NUMERICAL  TESTS  OF  THE  METHOD 


The  maximum  entropy  spectral  estimation  technique  was 
tested  by  representing  the  seismic  signal  as  a unit  spike  convolved  with  a 
luw-pass  filter.  Tests  were  performed  by  choosing  a corner  frequency  1.0  Hz 
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and  roil-off  of  12  dB /octavo.  The  system  response  was  approximated  by 
diRital  filters.  The  computed  spectrum  was  compared  to  the  theoretically 
expected  spectrum,  taking  into  account  effects  of  instrument  response  and 
filtering.  The  observed  spec  trum  was  constant  at  the  expected  level. 

A second  test  was  performed  to  investigate  the  influence  of 
the  coda  on  the  accuracy  of  spectra  computed  in  this  way.  The  test  data 
were  constructed  by  convolving  the  above  I Hz  corner  frequency  signal  gen- 
erated by  multiplying  a random  normal  variate  with  mean  zero  and  standard 
deviation  one  with  a factor  exp  (-N/50),  where  N was  the  index  of  the  point. 
Spectra  for  synthetic  signals  with  various  values  of  the  ratio  of  signal  energy 
to  coda  energy  (SCR)  as  determined  from  windows  picked  by  the  automatic 
detector,  arc  shown  in  Figure  III  - 1 . The  spectra  in  Figure  III  - 1 are  divided 
by  the  exact  theoretical  solution  for  the  signal  convolved  with  a unit  spike. 

Signal  to  coda  ratio  was  lowered  by  decreasing  the  spike  amp- 
litude. I he  calculated  spectral  amplitude  decreases  as  it  should.  Due  to  the 
effects  of  the  automatic  detector,  tne  expected  spectral  level  for  these  tests 
are  not  exactly  known.  Consequently  the  absolute  values  of  the  spectra  are 
not  significant.  The  relative  displacement  amplitude  levels  do  appear  to 
approximately  reflect  the  change  in  signal  amplitude.  Note  that  the  shape  of 
the  spectra  (divided  b/  the  expected  theoretical  solution  for  a spike)  should 
be  flat. 

For  SCR  - 4.0  the  spectrum  is  nearly  flat  as  expected,  ex- 
cept for  some  distortion  near  0.  7 Hz.  This  distortion  is  more  apparent  at 
SC  R 1.0.  At  SCR  0.  25,  the  technique  finds  negative  spectral  energy  in 
the  bands  near  0.7  Hz.  This  is  indicated  by  a discontinuous  plot  of  the  spectra 
with  negative  power  shown  as  blank  spaces.  Note  the  large  apparent  errors 
in  the  spectra  at  low  frequencies,  especially  near  blank  spaces. 

Ihe  presence  uf  this  divergence  suggests  that  as  the  coda 
energy  dominates  the  pulse  energy  in  the  seismic  signal  the  spectral  estimates 
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^ CX'><'c'|-d  »°  unreliable  at  low  frequencies.  Unfortunately 

this  situation  is  sometimes  encountered  with  real  seismic  data,  and  probably 
depends  primarily  on  the  size  of  the  coda  rather  than  event.  As  a result 
of  these  tests  i,  is  unknown  whether  nr  no.  our  spectra  are  reliable  at  low 
1 rcquencies.  If  reliable  spectra  are  to  be  obtained  at  all  frequencies,  or 
at  least  within  the  widest  possible  frequency  band,  errors  due  to  this  effect 
must  be  reduced  to  a tolerable  level. 


i 


!J 
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SECTION  IV 

RESULTS  WITH  SEISMIC  DATA 


Twenty  -eight  events  recorded  at  NORSAR  are  analyzed.  The 
event  parameters  are  listed  in  tW  IV-1.  All  magnitudes  are  determined 
teleseismicly  by  NOAA-PDE.  The  data  includes  6 regional  earthquakes  be- 
tween magnitude  3.1  and  4.5;  4 teleseismic  earthquakes  between  magnitude 
5.4  and  6.1;  2 regional  presumed  explosions;  and  8 teleseismic  presumed 
explosions  of  magnitude  between  4.  4 and  5.  6.  All  these  signals  have  signal- 
to-noise  ratios  greater  than  2.0. 

On  Figure  IV-1  are  earth  displacement  spectra  of  4 teleseismic 
presumed  explosions.  Most  spectra  exhibit  a broad  peak  at  about  2.  5 Hz  and 
fall  off  sharply  above  corner  frequencies  between  3.  2 and  4.  8 Hz.  At  lower 
frequencies  the  spectra  decrease  in  amplitude  at  a rate  of  0 to  6 dB/octave. 

The  effects  of  errors  in  the  displacement  spectrum  is  apparent  as  gaps  :n  the 
spectra.  These  occur  at  frequencies  where  negative  power  was  erroneously 
indicated.  The  seismometer  recordings  filtered  to  remove  poles  in  the  system 
response  are  shown  with  the  spectra. 

On  Figure  IV-2  are  earth  displacement  spectra  of  4 of  the 
regional  earthquakes.  Corner  frequencies  lie  between  2.  2 and  4.  5 Hz. 

Nearly  constant  amplitudes  appear  below  the  apparent  corner  frequencies. 
Considerable  effects  due  to  errors  in  calculating  the  ground  displacement 
spectra  are  indicated  by  gaps  in  tha  spectra,  and  by  large  and  erroneous 
variations  near  those  gaps.  Data  in  the  signal  window  are  also  shown. 

Where  the  signal-to-noise  ratio  is  large  it  is  possible  to 
compare  the  maximum  likelihood  spectra  with  conventional  spectra  computed 
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with  100  lags  of  the  autocorrelation  function  and  a Hanning  window.  Figure 
IV  -3  shows  such  a comparison  for  regional  earthquake  REQ3,  Figure  IV-4 
for  REQ12,  Figure  IV-5  for  REQ11;  Figure  IV-6  for  the  teleseismic  presumed 
explosion  TPE4;  Figure  IV-7  for  TPE3;  Figure  IV -8  for  the  teleseismic 
earthquake  TEQ6;  and  Figure  IV-9  for  TEQ4.  In  most  cases  the  corner  is 
indicated  with  less  ar  .biguity  by  the  new  method.  However,  t’is  se  are 
preliminary  results  and  should  be  treated  with  caution  until  the  effects  of 
errors,  whose  presence  is  obvious  from  inspection  of  the  results,  have  been 
resolved. 

Figure  IV-10  illustrates  three  regional  earthquakes  from  near- 
ly the  same  location  in  Greece.  These  events  are  gradually  emergent,  and 
ampl'/  demonstrate  the  ability  of  the  automatic  detector  to  pick  the  signal 
window. 

Corner  frequencies,  frequency  exponents,  and  displacement 
amplitudes  measured  from  earthquake  spectra  are  shown  in  Table  IV-2. 

A test  on  the  energy  was  made  to  check  whether  the  measurements  are  repre- 
sentative of  Drune's  model. 

In  Brunc's  model,  7 lies  between  1.0  and  2.0,  and  y ranges 

o 

from  2.0  to  10.0,  depending  on  the  fractional  stress  drop.  The  frequency 
exponents  shown  in  Table  IV-Z  appear  to  obey  these  requirements.  Conse- 
quently the  corner  frequencies  there  were  use  * to  estimate  € shown  on  Table 
IV-2  through  : 


€ = 0.8  log  (fQ-fc)  (I V -1) 

which  is  an  approximation  to  Brune's  theoretical  results. 

Equation  (II-4)  was  used  to  calculate  an  equivalent  corner 
frequency  under  Brune's  model  by  integrating  the  square  of  the  spectral 

density  to  find  the  radiated  seismic  energy  E . The  calculated  values  of  the 
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corner  frequency,  f , found  from  the  energy,  are  shown  in  Table  IV-2  and 

Figure  IV -11.  The  correlation  coefficent  between  f and  f i s 0.  95,  and  the 

e c 

standard  deviation  between  l^and  f is  0.083  H z.  Consequently  we  find  that 
there  is  good  consistency  between  the  model  and  the  data,  and  that  quantities 
other  than  corner  frequency  and  low  frequency  displacement  amplitude  can 
be  described  by  the  model. 

Spectral  magnitudes  estimated  by  an  equation  analogous  to 
equation  (II - 5)  for  local  magnitude,  are  given  by  equation  (IV-2). 


msb  = log  fio  + l"  5 lo&  fc  + k loR  A + C.  (IV-2) 

In  this  equation  C determined  empirically  for  the  earthquake  sample  by 

minimizing  the  variance  between  the  calculated  m , and  m measured  bv  the 

sb  b 7 

network.  The  factor  k for  distance  scaling  was  taken  to  be  + 3.  0 for  regional 
events  and  +1.0  for  teleseismic  events,  reflecting  the  presence  of  head  waves 
and  direct  waves  respectively.  For  the  regional  events,  an  average  value  of 
C was  found  to  be  -1.  39,  and  the  standard  deviation  between  the  magnitudes 
determined  by  equation  (IV-2)  using  this  value  and  the  network  magnitudes  was 
0.  17  magnitude  units.  For  the  teleseismic  events,  C was  +1.89  and  the 
standard  deviation  between  network  and  spectral  magnitudes  was  0.29  magni- 
tude units.  For  the  whole  sample,  the  standard  deviation  between  the  two 
magnitude  estimates  was  0.  20  magnitude  units,  and  the  two  magnitude  estimates 
were  highly  correlated  with  coeffecient  0.97. 

An  attenuation  factor 

C + k log  A = m - log  11  - — log  f 
b ° o 2 b c 

is  plotted  as  a function  of  epiccntral  distance  in  Figure  IV-12  for  each  event 
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EMPIRICALLY  DERIVED  SPECTRAL  AMPLITUDE 
SCALING  WITH  EPICENTRAL  DISTANCE 


in  the  earthquake  sample.  The  solid  lines  are  the  best  fit  to  the  data,  using 
the  values  of  C found  above.  The  transisition  from  regional  to  teleseismic 
behavior  occurs  at  39.  8°  for  this  data  set,  but  only  small  variations  in  the 
value  of  C would  allow  this  point  to  be  shifted  to  25°. 

To  investigate  its  power  as  a discriminant,  the  specral  mag- 
nitude computed  from  equation  (IV -2)  was  plotted  for  earthquakes  and  presumed 
explosions.  The  data  are  shown  on  Table  IV-2  and  plotted  on  Figure  IV-13. 

The  average  deviation  of  the  explosion  spectral  magnitudes 
from  their  network  magnitudes  was  0.  54  magnitude  units,  and  their  correla- 
tion with  the  network  m,  was  0.  98.  This  difference  is  reflected  in  the  dis- 

b 

placement  of  the  lines  representing  the  best  linear  fit  to  the  earthquakes  and 
presumed  explosions. 

Assuming  that  the  explosion  magnitudes  arc  normally  distri- 
buted about  their  regression  line,  95%  of  them  should  fall  above  the  dashed 
line.  All  the  earthquakes  fall  on  or  below  this  line,  and  only  one  presumed 
explosion  lies  below  it.  Furthermore,  that  explosion  is  believed  to  be  a 
multiple  event.  For  these  reasons,  this  criterion  has  promise  as  a discrimi- 
nant. It  should  be  emphasized,  however,  that  the  limited  data  sample  prevents 
hard  conclusions  from  being  drawn  on  this  subject. 

On  Figure  IV-14  is  a plot  of  measured  corner  frequencies  versus 
network  magnitude.  The  best  linear  fit  to  the  earthquake  population  is 


Log  f = - 0.074  m,  +0.74.  (IV-3) 

c b 

Six  of  the  logarithms  of  the  explosion  corner  frequency  measurements  lie 
further  than  0,  13  from  the  regression  line  for  Log  f . The  event  RPE5, 
which  failed  the  spectral  magnitude  versus  network  magnitude  discrimination 
test  was  claimed  to  be  an  explosion  by  this  corner  frequency  versus  network 
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FIGURE  IV  - 1 3 


MAGNITUDES  COMPUTED  FROM  SPECTRA  OF  NORSAR  DATA 
m.  COMPARED  WITH  EVENT  MAGNITUDES  AVERAGED 
FROM  TELESEISMIC  NOAA-PDE  MEASUREMENTS  OF  m 
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magnitude  test.  The  large  amount  of  scatter  evident  is  consistent  with  the 
results  of  Peppin  (1967)  and  Thatcher  and  Hanks  (1973). 


Figure  IV- 1 5 is  an  empirical  plot  of  Log  ft  corrected  for 
j . o 

distance  attenuation  defined  to  be 


Log  ftQ  = Log  ftQ  + k Log  \ + C 


(IV -4) 


where  C = -1.  39  and  k = 3 for  regional  events,  and  C = 1.  89  and  k = 1 for 

teleseismic  events.  The  median  line  for  the  average  Log  "o'  versus  net- 
work magnitude  is 


Log  ft  = 1.  14  m -1.21  . 
o b 


( I V-5) 


Another  linear  combination  of  Log^  and  Log  f , which  could 
provide  discrimination  information  is  the  logarithm  of  the  source  dimension, 
which  should  be  smaller  for  explosions  than  for  earthquakes.  From  equation 
(II-5),  Ix>g  (2r)  is  plotted  against  Am  = mbs-mb,  which  reflects  anomalous 
radiated  seismic  energy.  The  results  are  shown  on  Figure  IV-16,  where  all  of 
the  explosions  are  separated  except  RPE5.  Eight  of  the  ten  explosions  show 
well  defined  separation  from  the  space  occupied  by  earthquakes. 

Where  Alog  is  the  deviation  of  the  logarithmic  frequency 
from  its  expected  value  given  by  equation  (1V-3),  Figure  IV- 17  shows  Alog  f 
plotted  against  n^.  This  discriminant  combines  those  presented  in  Figures 
IV-13  and  IV-14.  Complete  discrimination  between  earthquakes  and  explosions 
is  obtained  by  the  solid  lines,  which  require  that  events  whose  corner  frequen- 
cy is  anomolously  high  or  whose  spectral  magnitude  is  anamolously  high  be 
classified  as  explosions.  However,  the  discriminant  based  solely  on  corner 

frequency  is  supported  by  only  one  point,  and  further  data  are  required  to 
establish  its  validity. 
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These  results  suggest  that  short-period  spectral  discriminants 
may  provide  at  least  supporting  evidence  for  classification  of  everts.  More 
development  of  the  maximum  entropy  spectrum,  and  a wider  data  base  might 
well  provide  a discriminant  powerful  enough  to  be  independent  of  other  tests. 
We  strongly  recommend  that  such  deve'opment  be  carried  out. 

Another  important  result  of  this  work  is  the  close  corres- 
pondence between  the  spectral  magnitudes  found  for  regional  events  and  the 
teleseismiciy  determined  magnitudes  (average  of  NOAA-PDE  magnitudes) 
of  these  same  events.  This  correspondence  suggests  that  the  spectral 
method  may  be  a convienent  and  accurate  means  of  determining  the  true  mag- 
nitude of  an  event  at  regional  distances.  The  absence  of  such  a simple  method 
for  regional  magnitude  determination  is  the  cause  of  a number  of  problems 
in  seismology  at  present. 
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SECTION  V 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  purpose  of  this  study  was  to  examine  sho  *t-period  spectral 
discriminants,  and  in  particular  discriminants  based  on  corner  frequency, 
which  is  expected  to  reflect  differences  in  source  dimensions.  In  order  to  do 
this  a method  of  estimation  based  on  th<i  maximum  entropy  spectrum  was  de- 
veloped. This  method  does  away  with  the  conventional  assumption  that  the 
data  are  zero  outside  the  sample  interval,  which  is  especially  unrealistic  for 
low  signal- to -noise  events.  This  result  is  achieved  by  estimating  the  auto- 
correlation function  outside  the  data  interval  with  a prediction  error  filter, 
and  deriving  the  spectrum  from  this  extended  autocorrelation  function.  Tests 
on  synthetic  data  showed  that  the  shapes  of  spectra,  and  possibly  the  ampli- 
tudes as  well,  are  quite  satisfactory. 

For  high  signal-to-noise  ratio  the  spectra  estimated  in  this  way 
are  in  general  agreement  with  those  calculated  by  the  conventional  Fourier 
transform  method,  but  show  less  ambiguity  in  their  corner  frequency.  At  low 
signal-to-noise  ratio,  Fourier  methods  cannot  estimate  these  parameters  at 

all.  while  the  maximum  entropy  method  yields  clear-cut  solutions  in  many 
cases. 

In  order  to  investigate  the  utility  of  the  method  as  applied  to  the 
discrimination  problem,  P-wave  maximum  e.  tropy  spectra  were  watched  for 
28  earthquakes  and  presumed  explosions. 

To  determine  the  internal  consistency  of  the  data,  and  to  find 
if  the  data  fit  the  model,  corner  frequencies  calculated  by  energy  considera- 
tions were  compared  with  those  measured  directly  from  spectra.  The  agree- 


ment was  quite  satisfactory. 


Corner  frequencies  measured  directly  from  spectra  were  plot- 
ted against  magnitude  for  the  whole  data  sample.  There  was  a general  trend 
for  events  with  larger  magnitudes  tc  have  lower  corner  frequencies,  but  the 
scatter  in  the  data  was  too  great  to  effectively  separate  earthquakes  and  pre- 
sumed explosions.  Consequently  we  conclude  that  corner  frequency  by  itself 
has  very  little  discriminatory  power. 

A spectral  magnitude,  whose  definition  was  motivated  by  that 
of  spectral  magnit  ides  calculated  for  local  events,  was  found  for  the  earth- 
quake sample.  It  was  a function  of  the  zero-frequency  spectral  amplitude, 
the  corner  frequency,  the  epicentral  distance,  and  one  constant  which  was  ad- 
justed to  minimize  the  variance  between  the  spectral  and  network  magnitudes. 
The  network  and  spectral  magnitudes  found  in  this  way  were  linearly  related 
with  a slope  of  nearly  one.  Then  spectral  magnitudes  were  found  for  the  pre- 
sumed explosions,  using  the  relationship  derived  for  earthquakes.  These 
magnitudes  were  significantly  higher  than  the  corresponding  network  magni- 
tudes, reflecting  the  higher  radiated  seismic  energy  associated  with  explo- 
sions. A discriminant  based  on  this  magnitude  difference  would  have  separat- 
ed all  earthquakes  from  explosions,  and  would  have  misclassified  only  one 
explosion. 

The  success  of  this  single-valued  discriminant  lead  to  the 
search  for  a more  powerful  test  using  two  parameters.  A discriminant  which 
appears  to  be  effective  combines  the  difference  between  spectral  magnitude 
and  conventional  event  magnitude  with  deviations  of  the  corner  frequency  from 
the  value  obtained  from  a regression  of  corner  frequency  versus  magnitude. 
This  two  dimensional  plot  of  short-period  discriminants  indicated  improved 
capability,  but  the  data  sample  was  still  too  small  to  draw  reliable  conclusions. 
The  results  suggest  that  a measure  of  source  dimension  has  discriminating 
power  and  that  further  research  on  this  topic  is  justified. 
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p J oblems  remain  with  maximum  entropy  spectra,  how- 
ever. In  some  situations  additive  random  noise,  coda,  roundoff  errors,  and 
inherent  limitations  of  the  algorithm  result  in  spectral  power  estimates  less 
than  zero  at  some  frequencies.  This  is  the  most  serious  problem  with  the 
maximum  entropy  specti  um  at  this  time. 

In  view  of  the  improved  spectra  derived  for  events  with  low 
signal-to-noise  ratio,  and  the  more  clear-cut  measurements  available  fror  1 
spectra  with  high  signal-to-noise  ratio  as  compared  with  conventional  tech- 
niques, we  recommend  that  more  work  on  this  subject  be  carried  out.  In  par- 
ticular, a more  accurate  maximum  entropy  spectral  estimation  technique 
should  be  developed.  The  improved  method  should  be  used  to  obtain  accurate 
displacement  spectrum  of  events  for  frequencies  greater  than  0.25  Hz.  Th«i 
method  should  be  capable  of  minimizing  errors  due  to  ambient  noise  and  coda. 
Given  such  an  improved  method,  it  should  be  applied  to  a broader  data  base 
comprising  both  teleseismic  and  regional  events.  The  magnitude  base  lines 
of  both  types  of  events  should  be  extended  as  far  aj  possible  to  determine  any 
magnitude  limitations  on  the  application  of  the  discriminant. 
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